#----generating multivairate Normal

library(MASS)
library(irr)
library(Hmisc) # for the errorbar plot

PercentageExactMatch <-function(trueScore,fittedScore){
  nmatch = length(trueScore[trueScore == fittedScore])
  ntotal = length(trueScore)
  return(nmatch/ntotal*100) 
}

roundScore3<-function(escore){
  d1 = escore< 2
  d2 = escore>=2 & escore<4
  d3 = escore >= 4
  escore[d1] = 1
  escore[d2]=2
  escore[d3] = 3
  return(escore)
}

roundScore4<-function(escore){
  d1 = escore< 2
  d2 = escore>=2 & escore<3
  d3 = escore >= 3 & escore<4
  d4 = escore >= 4
  escore[d1] = 1
  escore[d2]=2
  escore[d3] = 3
  escore[d4] = 4
  return(escore)
}


roundScore5<-function(escore){
  d1 = escore< 1.5
  d2 = escore>=1.5 & escore<3
  d3 = escore >= 3 & escore < 4
  d4 = escore >= 4 & escore < 5.5
  d5 = escore>=5.5
  escore[d1] = 1
  escore[d2]=2
  escore[d3] = 3
  escore[d4] = 4
  escore[d5] = 5
  return(escore)
}


roundScore6<-function(escore){
  d1 = escore< 1.5
  d2 = escore>=1.5 & escore<2.5
  d3 = escore >= 2.5 & escore < 3.5
  d4 = escore >= 3.5 & escore < 4.5
  d5 = escore>= 4.5 & escore < 5.5
  d6 = escore>=5.5
  escore[d1] = 1
  escore[d2]=2
  escore[d3] = 3
  escore[d4] = 4
  escore[d5] = 5
  escore[d6] = 6
  return(escore)
}


roundScore7<-function(escore){
  d1 = escore< 1
  d2 = escore>=1 & escore<2
  d3 = escore >= 2 & escore < 3
  d4 = escore >= 3 & escore < 4
  d5 = escore>= 4 & escore < 5
  d6 = escore>=5.5 & escore < 6
  d7 = escore>=6
  escore[d1] = 1
  escore[d2]=2
  escore[d3] = 3
  escore[d4] = 4
  escore[d5] = 5
  escore[d6] = 6
  escore[d7] = 7
  return(escore)
}

exactMatchCorrelation<-function(correlation,nlevel){
  sigma = matrix(1,nrow=2,ncol=2)
  sigma[1,2] = correlation
  sigma[2,1] = correlation
  t=mvrnorm(n = 1000, mu=c(3,3), Sigma=sigma)
  h1 = t[,1]
  h2 = t[,2]
  if (nlevel == 3){
    h1 = roundScore3(h1)
    h2 = roundScore3(h2)
  }
  if (nlevel == 4){
    h1 = roundScore4(h1)
    h2 = roundScore4(h2)
  }
  if (nlevel == 5){
    h1 = roundScore5(h1)
    h2 = roundScore5(h2)
  }
  if (nlevel == 6){
    h1 = roundScore6(h1)
    h2 = roundScore6(h2)
  }
  if (nlevel == 7){
    h1 = roundScore7(h1)
    h2 = roundScore7(h2)
  }
  pmtch = PercentageExactMatch(h1,h2)
  correlationCal = cor(h1,h2)
  kapa = kappa2(cbind(h1,h2),weight="squared")
  return(c(correlationCal,pmtch,kapa$value))
}

generateResult<-function(nlevel=3){
  correlation = seq(0.1,1,0.05)
  n = length(correlation)
  m = 100
  exactMtch = matrix(0,nrow=n,ncol=m)
  corrNew = matrix(0,nrow=n,ncol=m)
  kapaNew = matrix(0,nrow=n,ncol=m)
  for(i in 1:n){
    for (j in 1:m){
      t = exactMatchCorrelation(correlation[i],nlevel)
      corrNew[i,j] = t[1]
      exactMtch[i,j] = t[2]
      kapaNew[i,j] = t[3]
    }
  }
  MatchPercentMean = rowMeans(exactMtch)
  MatchPercentSd = apply(exactMtch,1,sd)
  corrMean  = rowMeans(corrNew)
  corrSd = apply(corrNew,1,sd)
  kapaMean = rowMeans(kapaNew)
  kapaSd = apply(kapaNew,1,sd)
  
  data = cbind(corrMean,kapaMean,MatchPercentMean)
  write.table(data,file=paste('~/research/LPworknew/data_',as.character(nlevel),'_level.csv',sep=''),sep=',',row.names=F,col.names=T)
  
  png(paste('~/research/LPworknew/correlation_exactmatch_',as.character(nlevel),'level.png',sep=''))
  errbar(corrMean,MatchPercentMean,MatchPercentMean+MatchPercentSd,MatchPercentMean-MatchPercentSd,xlab='correlation',ylab='Percentage of exact match')
  text(0.2,90,paste('# of score levels:',as.character(nlevel)))
  grid()
  dev.off()
  
  png(paste('~/research/LPworknew/kappa_exactmatch_',as.character(nlevel),'level.png',sep=''))
  errbar(kapaMean,MatchPercentMean,MatchPercentMean+MatchPercentSd,MatchPercentMean-MatchPercentSd,xlab='quadratic weighted kappa',ylab='Percentage of exact match')
  text(0.2,90,paste('# of score levels:',as.character(nlevel)))
  grid()
  dev.off()
  
}

#---call functions-----

generateResult(3)
generateResult(4)
generateResult(5)
generateResult(6)
generateResult(7)



